Optimized random phase approximations for arbitrary reference systems: 
extremum conditions and thermodynamic consistence 
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The optimized random phase approximation (ORPA) for classical liquids is re-examined in the 
framework of the generating functional approach to the integral equations. We show that the two 
main variants of the approximation correspond to the addition of the same correction to two different 
first order approximations of the homogeneous liquid free energy. Furthermore, we show that it 
is possible to consistently use the ORPA with arbitrary reference systems described by continuous 
potentials and that the same approximation is equivalent to a particular extremum condition for the 
corresponding generating functional. Finally, it is possible to enforce the thermodynamic consistence 
between the thermal and the virial route to the equation of state by requiring the global extremum 
condition on the generating functional. 
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I. INTRODUCTION 

The optimized random phase approximation (ORPA) Q| for classical liquids has been extensively used in the last 
decades to obtain information on the structure and, to a minor extent, on the thermodynamics of simple liquids 
and mostly liquid metals The approximation was originally developed in the context of the perturbative 

approach to the thermodynamic and structure of simple liquids . More recently, this approximation has been used 
as ingredient for studying liquids in porous media and critical phenomena in simple liquids [^. The standard 
implementation of ORPA is based on the splitting of the interparticle potential into a repulsive (reference) and an 
attractive (perturbation) part. The effect of the attraction on the pair correlation functions of the purely repulsive 
reference system is treated at the level of the random phase approximation (RPA) at large distances and by enforcing 
the excluded volume effect at short distances. 

Numerical studies have shown that the ORPA yields a very accurate description of the structure factor and ther- 
modynamics of simple fluids. At present, for simple liquids, the quality of the ORPA results is comparable to that of 
. , statc-of-thc-art calculations based on the modified hypernetted chain (MHNC) approximation [^,0 or other modern 
integral equations like the HMS A . 
■ ■ Usually, the reference system is modeled by hard sphere interactions. However, in some cases, either the nature 

of the interactions or results from other theoretical approaches, indicate that a soft reference system could provide a 
better reference system for the ORPA. The original derivation of the ORPA does not allow a direct extension of the 
formulae to the case of a reference system interacting without hard core. A non-standard implementation for liquid 
alkali metals using a one component plasma as reference system Q, although providing good results, was not general 
enough and was subject to some criticisms JT2[ . 

More recently, motivated by the need of improving some variational calculations for liquid metals we reviewed 
the ORPA from the computational as well from the theoretical point of view. In a previous paper p4| , we showed 
that the solution of the ORPA equations is unique and we proposed a new robust, accurate and efficient numerical 
algorithm to solve the equations. In the present paper we address the problem of reformulating the theory in order 
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to deal in a consistent way with continuous reference systems without any intermediate introduction of auxiHary 
hard-sphere systems Q. We give a solution to this problem in the same spirit of Gillan's extension of the mean 
spherical approximation for soft potentials and subsequent elaborations by Rosenfeld [|l6| . 

We found useful to recast the ORPA equations in the framework of the generating functionals for the integral 
equations of the theory of liquids We show that the two different first-order expansions of the free energy 

functional corresponding to the well known Gibbs-Bogolioubov and Weeks-Chandler and Andersen approximations 
for the free energy |]l8| can be transformed into two closely related forms of the ORPA by addition of the same 
functional. The resulting correlation functions differ by the choice of the reference system pair correlation function. 

Moreover, we are able to show that our condition for a continuous ORPA correction to the pair correlations is 
equivalent to an extremum condition for the ORPA generating functional. Since only the variational determination 
of all the free parameters of the functional allows the identification of its value with the Helmoltz free energy, this 
choice enforces the thermodynamic consistence in a natural way. 

The paper is organized as follows. In section 2, we show how two versions of the ORPA differing only in the 
treatment of the reference system correlations can be obtained from a variational problem for two related functionals. 
In section 3 we show that it is possible to define a consistent ORPA for continuous potentials and that such an 
extension is equivalent to an additional requirement of extremum for the ORPA functionals. In section 4 the issue of 
the thermodynamic consistence of the two approximations is briefly discussed. Conclusions are summarized in section 
5. 



II. TWO GENERATING FUNCTIONALS FOR THE ORPA 



The starting point of the ORPA is a suitable decomposition of the interatomic potential (j>{r) into a reference 
potential 4>o{r) and a perturbation (the rest) 4>i{r): 

(t>{r) = Mr) + Mr) (1) 

Although the original ORPA was based on a specific choice of such a decomposition, in the following discus- 
sion we temporarily leave unspecified the exact characterization of (j^oir)- Equation (1) naturally leads to a similar 
decomposition of the total and the direct correlation functions h(r) and c(r): 

h{r) ^ ho{r) + Ah{r) (2) 
c(r) =co(r) + Ac(r) (3) 

where hQ(r) and co(r) are the correlation functions of a reference fluid whose particles interact via the potential 
4>o{r). The thermodynamics and the correlation functions of the reference system are considered as known quantities. 
Ah{r) and Ac(r) are deflned by equations (2) and (3) and arc the unknown functions of the theory. A relation among 
them, for a fluid whose number density is p, is provided by the Ornstein-Zernike equation 

h{r) = c(r) +p [ d^r'h{r')c{\r - r'|) (4) 



which, taking into account the fact that h(){r) and co(r) do satisfy the same equation, results in the following 
relation between the Fourier transforms of Ah{r) and Ac(r): 

1 - pAc{q)So{q) 

In formula (j^), So{q) = 1 + pho{q) is the structure factor of the reference system. A hat on a function of q indicates 
the three-dimensional Fourier transforms of the corresponding function defined in the r-space while p is the number 
density of the system. 

So far no approximation has been introduced yet. By complementing equation with any approximate relation 
between Ah{r), Ac{r) and 0i(r) we get a closed set of non-linear integral equations that has to be solved. 
In particular, the ORPA closure corresponds to the dual relations: 

Ac(r) = —(3(j)i{r) for r > cr (6) 
Ah{r) =0 for r < cr . (7) 
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These equations impose, up to the finite crossover distance a, the matching of the asymptotic long range behavior 



of Ac(r) (eq. (|6 
distances (eq. ( 



and the condition that the approximation would not modify the pair correlation function at short 
). Due to the presence of relation one could use as independent variable either the values of 
Ah{r) at distances beyond a or, more conveniently, the values of the function x(^) = Ac(7') for r < a (xif) = for 
r > 0) . 

It is easy to show that equation (0) corresponds to the extremum condition for the following functional []l|jl7| of 
X{r): 



FRjNcixir)] = 



2(27r)V 



d'q{log[l + So{q)p{q)]-p{q)So{q)} 



(8) 



where p{q) = p{(3(j>i{q) — x{q))- ^ diagrammatic treatment, Fjunq would correspond to the sum of rm^-like 
diagrams and, as it is well known, for x(r) = corresponds to the random phase approximation which usually 
violates the core condition (|^). The ORPA enforces such a condition. 

Indeed, by taking the functional derivative with respect to Ac{q) we have: 
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-Ah{q) 



SAc{q) (27r)= 

and Fourier-transforming to the r-space we get for all the values of r: 

SFrjng 



(SAc(r) 



Ah{r). 



(9) 
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For < r < a, equation (|l0| ) becomes an integral equation for the unknown %('')) different from zero only in such 



a region: 
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for r < a 



(11) 



Thus, equation (11) is equivalent to imposing an extremum condition on FRjf^fQ with respect to variations of 
Ac(r) — x(r) (inside a). If the reference potential is such that goir) inside a is zero, we see that the extremum 
condition is equivalent to the physical requirement that the size of the exclusion hole of the reference system is 
preserved by the perturbation. 

It is quite easy to verify that the solution of the equation (|^), provided it exists, is actually unique and corresponds 
to the maximum of the ORPA generating functional jl^ (||). 

To complete the description of the system, an explicit prescription for the reference system pair correlation function 
go{r) is required. In the usual approach to ORPA, the choice of gQ{r) is treated as a separate step. Here we prefer 
to define a "total" generating functional from which the full ORPA g{r) is derived. Actually we can introduce two 
functionals having both Frimg as generator of the "ORPA" contribution to the pair correlations and differing in the 
resulting g^. 

For a homogeneous liquid interacting through a pair potential (/>, the Helmoltz free energy per particle F can be 
considered [M a functional of 0(r) as well as a functional of the function e(r) = e~^'^'^^\ It is easy to show that 



5F 2 
)<j)[r) p 



and that 



5F 



Se{r) p 

where y{r) is the so-called cavity correlation function y(r) 
and J-^wca[g] — as follows. 



(12) 



(13) 



g{r)e'^'^^\ We introduce two functionals — TGB[A(j)] 



Tgb[A^] = ^ / drgo(r)/3A0(r 



-(- 

2p \2tt 



dq 



p{q)Soiq) ~ logil+p{q)So{q)) 



(14) 



^wcA [e] 



dryo(r)Ae(r) 
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2p \2tt 



dq 



p{q)So(q) - log(l+p{q)So{q) 



(15) 
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By functional differentiation of Tgb and TwcA with respect to 4>{r) and e(r) respectively, we get 



g = 9o + Ah{r) (16) 

and 

y^yo + e^-^'Ahir). (17) 
From the last equation we get immediately the g{r) resulting from J-wcA as 

9 = 9oe-^^^^^ + Ah{r). (18) 

Thus, the functionals J-qb and J-'wcA are such that the deviation from the reference system pair correlation function 
is always given by the ORPA approximation Ah{r) (eqn. but the reference system pair correlation function is 
9q in one case and yos~^'^^ in the other case. Due to the form of the reference system pair correlation functions 
and the corresponding generating functionals, we refer to the former approximation as the Gibbs-Bogoliubov ORPA 
(GB-ORPA) and to the latter as the Weeks-Chandler- Andersen ORPA (WCA-ORPA). 

Notice that at this level the two functionals have been introduced just as generating functionals for the pair 
correlation functions and we are not allowed yet to identify the values of the two functionals at the extremum with 
the Helmoltz free energy. 



III. ORPA FOR CONTINUOUS POTENTIALS 

For a general value of the parameter a, the solution x(r) of equation ( [III ) and the resulting Ah{r) are discontinuous 
at a no matter if the reference system potential is continuous or not. While such a discontinuity looks relatively 
harmless if the reference 9o{r) has a hard core of diameter cr, a discontinuity in Ah{r) at would be spurious in 
connection with a continuous reference system. 

For a similar problem, occurring in the case of the mean spherical approximation (MSA), a satisfactory solution 
was found [ p5[ by determining a in such a way that the resulting correlation functions were continuous at a. Here, 
we can similarly impose the continuity of Ac(r) (or equivalently Ah(r)) at r = cr. Thus , we add the condition 

Ah{<T+) ^ (19) 

as additional equation for a. 

Moreover, still in analogy with the MSA case, we can prove that the continuity condition at a is equivalent to an 
extremum condition of the ORPA functional (^ as a function of a. 
As shown in appendix A we have 

dFoRPA/da = 27:a^Ax^{a-). (20) 

Thus, the continuity condition on the correlation functions implies that the GB-ORPA and WCA-ORPA functionals 
have an extremum (inflection point) at a. As we will discuss in the next section, this extremum condition is also the 
clue for a thermodynamic consistent theory. 

Here we just notice that there is a manyfold of solutions of equation (|l^). However a lower limit for a is given by 
the size of the excluded volume region of the reference system. That is, the region such that 

9{r) « (21) 

A choice of a smaller than the reference system exclusion hole would result again in an unphysical ORPA g{r). On 
the other hand, since ^wCA and J-^gb are increasing functions of a (eqn. |2^), the minimum value of will be achieved 
for the first value of a larger than the reference system exclusion hole. Moreover, increasing a, the size of the ORPA 
correction to the reference system thermodynamics and correlations rapidly decreases. 



IV. GENERATING FUNCTIONALS AND THERMODYNAMIC CONSISTENCE 



Now we are in the position to discuss the thermodynamic interpretation of the functionals J-qb and J-'wcA and 
the specific issue of the thermodynamic consistency. 
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It is well known that approximate integral equation theories for the correlation functions show quantitative violations 
of fundamental thermodynamic equalities. In particular, here we are concerned with the equalities generated by the 
identification of the generating functional with the Helmoltz free energy per particle /. The most obvious of such 
equalities is the equality between the pressure p obtained from the free energy per particle /, 



pp dm 

— = P 22) 

P P 



and that found through the virial theorem, 



— [ 9{r)r(3q^'{r)dr (23) 

p 6 J 

A necessary condition to ensure that a functional F[<j)] is actually a free energy functional is the validity of equation 
(H) (or (©) @. ^ ^ 

Such a condition would be fulfilled by the functionals defined in equations lR4) and ( p5| ) if the dependence of such 
functionals on all the parameters of the reference system, on a (say a^) and on vanishes. Then we have to satisfy 
the following equations: 



dF 

dai 

dF 
_5F_ 



^ (24) 
= (25) 

(26) 



Eqn (26) corresponds to the ORPA formula [7] while eqn (24) is a way of determining the reference system param- 
eters. Eqn. (|2^), as we have shown in the previous section is also related to the continuity of the resulting correlation 
functions. Therefore, in order to have thermodynamic consistency we have to ensure that the functional would be 
extremum with respect to variations of all the parameters. An analogous requirement for the choice of the reference 
system in connection with the modified hypernetted chain approximation (MHNC) was derived by Lado et al. 
Even closer to the present problem is the analysis of the choice of the reference system within the WCA perturbation 
theory provided by Lado . 

Different choices of the reference system parameters are certainly conceivable and actually this is the existing 
situation. It is not easy to anticipate what is the best choice for all possible systems and a final assessment should be 
left to explicit numerical investigations. However, here we can notice that only the choices corresponding to extrema 
of the generating functionals or choices completely independent of the thermodynamic state would ensure the free 
energy nature of the generating functionals and then, as a consequence, the consistence of the energy and virial routes 
to the equation of state. 

V. CONCLUSIONS 

In the present paper we have rephrased the ORPA in the language of the generating functionals for the pair 
correlation function. In this way we could easily obtain three main results: 

1) we can derive from a unified treatment the two prescription for the reference system g{r) ( equations ( |l6| ) and 
( p^ ) present in the literature; 

2) we can show how the ORPA can be extended to deal with continuous reference system interactions, potentially 
increasing the range of applicability of this approximation; 

3) we show that the closure equations, the removal of the discontinuity in the resulting pair correlations and the 
identification of the generating functionals with the Helmoltz free energy can be reduced to the unique and 
unifying requirement of a variational principle on the functionals with respect to all the independent variables 
and parameters. 
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The theory presented in this paper provides a general scheme corresponding to many possible choices for the 
individual ingredients of the ORPA. Actually, depending on the reference system and on the flavor of the ORPA (GB 
or WC A) , we have introduced different possibilities. For this reason we postpone detailed numerical investigations to 
the application of the approximation to specific problems. 

Taking into account the already satisfactory level of accuracy of the standard implementations of the ORPA, and 
judging from preliminar calculations, we can anticipate a good quality of the resulting numerical results. In particular 
thermodynamical investigations could now benefit from the clarified status of thermodynamic consistency in the 
ORPA. In this respect, we believe that the ORPA could play an important role as one of the best candidate for the 
investigation of the fluid phase diagrams. 
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APPENDIX A 

The proof of eqn. (^) is given as follows: 
Let 

p{q) PpMq) - px{q)- (Al) 

Since the only dependence of J-wcA or J-wcA on a is through x, the derivative of Fujng (eqn. (^), we have to 
evaluate 



dFRjNG 1/1 
da ~ 2p V27r 



dq 



Soiq) - S{q) 



Now, taking into account the finite support of x(^)) 



— — cx(c )stn[qa) + / r — sin(qr)dr. 



da q Jf) da q 

By using Parseval's equality, eqn. (A3) becomes: 



OFring 1 f ^ Mr) If 1 



3 /.oc 

2y d^^M^)^^!^ - ^ ) (47r)2y^ qA/i(g)CTx('^)sm(ga)dg (A4) 

the first term in eqn. (A4) is zero because when Ah ^ 0, the other term is zero and the reverse also follows. Eqn 
(A4) eventually reduces to 

^^1^;^ = -2na\ia-)Ahia+) = -2^a\^a-) (A5) 



giving eqn. (Ed)) when we take into account that Fmf^G appears in eqns. (|14|) and (15) with a negative sign 
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